#ifndef _CFDTD3D_
#define _CFDTD3D_

#define MaxMedium 500
#define Eps0 8.854E-12
#define Miu0 1.256637E-6
#define Pi 3.1415926
#define C0 3.0E8
#define Boltzmann_k 1.38E-23

//int main(int grac, char *grav[]);
int ReadData_Init();
int MainFDTD();
void CalculateD();
void AddSourceD();
void BoundD_2_Zero();
void AddSource();
void AddSourceE();
void SavePortData();
void SaveFieldonPlane();
void CalculateE();
void ConductorE_2_Zero();
void BoundE_2_Zero();
void CalculateH();
void AddSourceH();
void BoundH_2_Zero();
void DelData();
void PostProcessing();
void PostCalculateZ();
void PostCalculateS();
void PostCalculateBeita();
void AddSoftSource();

//Field
int MaxX=100,MaxY=100,MaxZ=100;
double dx,dy,dz,dt,dx_rate,dy_rate,dz_rate,min_dxyz;
double ***Dx, ***Dy, ***Dz;
double ***Ex, ***Ey, ***Ez;
double ExMax, EyMax, EzMax, EzMax1, EzMax2;
double ***gax, ***gay, ***gaz;
double ***gbx, ***gby, ***gbz;
double ***gcx, ***gcy, ***gcz;
double ***gdx, ***gdy, ***gdz;
double ***Ix, ***Iy, ***Iz;
double ***Sx, ***Sy, ***Sz;
double ***Hx, ***Hy, ***Hz;
double ***Idxl, ***Idyl, ***Idzl;
double ***Idxh, ***Idyh, ***Idzh;
double ***Ihxl, ***Ihyl, ***Ihzl;
double ***Ihxh, ***Ihyh, ***Ihzh;
double ***HB;
double Dtemp, curl_h, curl_e;

//Medium Para
double Epsilon_r, Epsilon_s, Miu_r, Sigma, DebyeTao;
int Conform, ConType;
int ***Model;
double Medium_EpsilonR[MaxMedium],Medium_EpsilonS[MaxMedium];
double Medium_MiuR[MaxMedium],Medium_Sigma[MaxMedium],Medium_DebyeTao[MaxMedium];
int Medium_Conform[MaxMedium],Medium_ConType[MaxMedium];
int NumMedium, MediumID;

//Conduction Plane
int NumNormalX, NumNormalY, NumNormalZ;
int **NormalX, **NormalY, **NormalZ;

//PML
double *gx3, *gx2, *gx1, *fx1, *fx2, *fx3;
double *gy3, *gy2, *gy1, *fy1, *fy2, *fy3;
double *gz3, *gz2, *gz1, *fz1, *fz2, *fz3;
int PML_lx, PML_ly, PML_lz;
double xnum,xd,xxn,xnpml,ynum,yd,yyn,ynpml,znum,zd,zzn,znpml;
double PML_factor=0.333,PML_m=4.,PML_epsr=8.;

//Source
int SourceType,SourceHS,T;
int SourceX1,SourceY1,SourceZ1,SourceX2,SourceY2,SourceZ2;
int SourceNormal,SourceVector;
double UWB1K=1, UWB1E0=1;
double UWB1tao,UWB1delta,UWB1t0,UWB1tr,UWB1tp;
double UWB2K=1;
double UWB2alpha,UWB2beita,UWB2t0,UWB2E0;
double Gausst0=0, GaussSpread=0;
double SinFreq=0;
double SourcePluse;
double RecSint1,RecSint2,RecSinf0;
//Plane Wave Source
double *E_inc,*H_inc,Epsr_Source;

//Watch
int NumWatch;
int *WatchX, *WatchY, *WatchZ;

//Port to save
int NumPort,iPort;
int *PortVX1,*PortVY1,*PortVZ1,*PortVX2,*PortVY2,*PortVZ2;
int *PortIX1,*PortIY1,*PortIZ1,*PortIX2,*PortIY2,*PortIZ2;
double Vtemp,Itemp1,Itemp2;
int PortVector,PortNormal,ixtemp,jytemp,kztemp;

//Field on Plane to Save
int NumFieldPlane,iFP;
int *FieldPlaneNormal,*FieldPlaneVector,*FieldPlaneTimeStep;

//Save Data
FILE *fInput, *fModel;
FILE *fWatch;
FILE *fSource,*fVPort,*fIPort1,*fIPort2;
FILE **fFieldPlane;

//Flag
int isRunning=0;
int OpenFileFlag=0,ReadDataFlag=0,DelDataFlag=0;
//int ExitFlag=0;
//int RefreshStep=1;

int StepStart, nsteps;
int i,j,k,n;
time_t t_start,t_end;

/*//PostProcessing
double fstart,fend,freq,df;
double *Z_r,*Z_i,**S_r,**S_i,*Beita_r,*Beita_i;
double **Vt_r,**Vt_i,**It1_r,**It1_i,**It2_r,**It2_i;
double **Vw_r,**Vw_i,**Iw1_r,**Iw1_i,**Iw2_r,**Iw2_i;*/

#endif
